///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file formic/vector/matrix.h
///
/// \brief   header file for the formic::matrix class
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef FORMIC_VECTOR_MATRIX_HEADER
#define FORMIC_VECTOR_MATRIX_HEADER

#include <iostream>
#include <string>

#include <boost/format.hpp>

#include <formic/vector/vector.h>
#include <formic/numeric/numeric.h>
#include <formic/lapack/interface.h>

namespace formic {

  // pre-declaration of non-const matrix
  template<class S, class L> class Matrix;

  // S is data type, L is length type
  template<class S, class L> class ConstMatrix {

    protected:

      L _n;

      L _m;

      S * _p;

      boost::shared_ptr<formic::VectorLock<S,L> > _lock;

    protected:

      void initialize_with_new_lock(const L n, const L m) {
        if (n < 0)
          throw formic::Exception("cannot initialize a matrix with a negative number of rows");
        if (m < 0)
          throw formic::Exception("cannot initialize a matrix with a negative number of columns");
        _n = n;
        _m = m;
        _p = 0;
        _lock.reset();
        if ( n*m > 0 ) {
          boost::shared_ptr<formic::VectorLock<S,L> > new_lock( new formic::VectorLock<S,L>(n*m) );
          _lock = new_lock;
          _p = _lock->cptr();
        }
      }

      void initialize_without_lock(const L n, const L m, const S * const p) {
        if (n < 0)
          throw formic::Exception("cannot initialize a matrix with a negative number of rows");
        if (m < 0)
          throw formic::Exception("cannot initialize a matrix with a negative number of columns");
        if ( n*m > 0 && p == 0 )
          throw formic::Exception("cannot initialize a matrix of nonzero size with a null pointer");
        if ( n*m == 0 && p != 0 )
          throw formic::Exception("cannot initialize a matrix of zero size with a non-null pointer");
        _n = n;
        _m = m;
        _p = (S *)p;
        _lock.reset();
      }

    private:

      // assignment operator disabled for constant matrix
      ConstMatrix & operator=(const ConstMatrix & other);

    //////////   constructors   //////////
    public:

      // constructor without initialization of elements (also the defualt constructor)
      explicit ConstMatrix(const L n = 0, const L m = 0) {
        this->initialize_with_new_lock(n, m);
      }

      // constructor with initialization of elements
      explicit ConstMatrix(const L n, const L m, const S v) {
        this->initialize_with_new_lock(n, m);
        std::fill(_p, _p + _n*_m, v);
      }

      // constructor with external storage and without initialization of elements
      explicit ConstMatrix(const L n, const L m, const S * const p) {
        this->initialize_without_lock(n, m, p);
      }

      // copy constructer (shallow copy)
      ConstMatrix(const ConstMatrix & other) {
        _n = other._n;
        _m = other._m;
        _p = other._p;
        _lock = other._lock;
      }

    //////////   element access   //////////
    public:

      const S & at(const L i, const L j) const {
        if ( i < 0 || i >= _n ) throw formic::Exception("out-of-range row index i = %i in formic::ConstMatrix::at for matrix of dimension %i by %i") % i % _n % _m;
        if ( j < 0 || j >= _m ) throw formic::Exception("out-of-range col index j = %i in formic::ConstMatrix::at for matrix of dimension %i by %i") % j % _n % _m;
        return _p[i+_n*j];
      }

      const S & operator()(const L i, const L j) const {
        #ifndef NDEBUG
        if ( i < 0 || i >= _n ) throw formic::Exception("out-of-range row index i = %i in formic::ConstMatrix::operator() for matrix of dimension %i by %i") % i % _n % _m;
        if ( j < 0 || j >= _m ) throw formic::Exception("out-of-range col index j = %i in formic::ConstMatrix::operator() for matrix of dimension %i by %i") % j % _n % _m;
        #endif
        return _p[i+_n*j];
      }

    //////////   other public members   //////////
    public:

      L size() const { return _n * _m; }

      L rows() const { return _n; }

      L cols() const { return _m; }

      std::string print(const std::string & fmt, const std::string & name = "") const;

      formic::Matrix<S,L> clone() const;

      formic::Matrix<S,L> operator-() const;

      formic::Matrix<S,L> operator+() const;

      formic::Matrix<S,L> t() const;

      formic::Matrix<S,L> c() const;

      template<class V>
      formic::Matrix<S,L> slice(const V & row_indices, const V & col_indices) const;

      // deep copy to c-style array
      const ConstMatrix & operator>>=(S * const x) const {
        std::copy(_p, _p + this->size(), x);
        return *this;
      }

      // deep copy to std::vector
      const ConstMatrix & operator>>=(std::vector<S> & x) const {
        if ( this->size() != x.size() )
          x.resize(this->size());
        if ( this->size() > 0 )
          *this >>= &x.at(0);
        return *this;
      }

      const S * begin() const { return _p; }
      const S * end()   const { return _p + this->size(); }

  };




  template<class S, class L> class Matrix : public formic::ConstMatrix<S,L> {

    public:

      // constructor without initialization of elements
      explicit Matrix(const L n = 0, const L m = 0) : formic::ConstMatrix<S,L>(n, m) {}

      // constructor with initialization of elements
      explicit Matrix(const L n, const L m, const S v) : formic::ConstMatrix<S,L>(n, m, v) {}

      // constructor with external storage and without initialization of elements
      explicit Matrix(const L n, const L m, S * const p) : formic::ConstMatrix<S,L>(n, m, p) {}

      // constructor with external storage and with initialization of elements
      explicit Matrix(const L n, const L m, S * const p, const S v) : formic::ConstMatrix<S,L>(n, m, p) {
        (*this) = v;
      }

      // copy constructer (shallow copy)
      Matrix(const Matrix & other) {
        this->_n = other._n;
        this->_m = other._m;
        this->_p = other._p;
        this->_lock = other._lock;
      }

      // assignment operator (shallow copy)
      Matrix & operator=(const Matrix & other) {
        this->_n = other._n;
        this->_m = other._m;
        this->_p = other._p;
        this->_lock = other._lock;
        return *this;
      }

      // deep copy to c-style array
      const Matrix & operator>>=(S * const x) const {
        std::copy(this->_p, this->_p + this->size(), x);
        return *this;
      }

      // deep copy from c-style array
      Matrix & operator<<=(const S * const x) {
        std::copy(x, x + this->size(), this->_p);
        return *this;
      }

      // deep copy to std::vector
      const Matrix & operator>>=(std::vector<S> & x) const {
        if ( this->size() != x.size() )
          x.resize(this->size());
        if ( this->size() > 0 )
          *this >>= &x.at(0);
        return *this;
      }

      // deep copy from std::vector
      Matrix & operator<<=(const std::vector<S> & x) {
        if ( this->size() != x.size() )
          throw formic::Exception("illegal attempt to deep copy from a std::vector to a formic::Matrix whose sizes do not match");
        if ( x.size() > 0 )
          *this <<= &x.at(0);
        return *this;
      }

      // deep copy from another formic::Matrix
      Matrix & operator<<=(const formic::ConstMatrix<S,L> & other) {
        if ( this->size() == other.size() ) {
          this->_n = other.rows();
          this->_m = other.cols();
        } else {
          if ( this->rows() * this->cols() > 0 && ! (this->_lock) )
            throw formic::Exception("illegal attempt to deep copy to a matrix with external storage that is not the correct size");
          this->initialize_with_new_lock(other.rows(), other.cols());
        }
        if ( this->size() > 0 )
          *this <<= &other.at(0,0);
        return *this;
      }

      /// \brief  clears the matrix, leaving it as a 0 by 0 empty matrix
      void clear() { *this = formic::Matrix<S,L>(); }

      /// \brief  clears the matrix and re-allocates it to be the specified size
      Matrix & reset(const int n, const int m) {
        this->clear();
        return *this = formic::Matrix<S,L>(n,m);
      }

      /// \brief  clears the matrix and re-allocates it to be the specified size, initializing its elements to v
      Matrix & reset(const int n, const int m, const S & v) {
        this->clear();
        return *this = formic::Matrix<S,L>(n,m,v);
      }

      const S & at(const L i, const L j) const {
        const int n = this->_n;
        const int m = this->_m;
        if ( i < 0 || i >= n ) throw formic::Exception("out-of-range row index i = %i in formic::Matrix::at for matrix of dimension %i by %i") % i % n % m;
        if ( j < 0 || j >= m ) throw formic::Exception("out-of-range col index j = %i in formic::Matrix::at for matrix of dimension %i by %i") % j % n % m;
        return this->_p[i+n*j];
      }

      S & at(const L i, const L j) {
        const int n = this->_n;
        const int m = this->_m;
        if ( i < 0 || i >= n ) throw formic::Exception("out-of-range row index i = %i in formic::Matrix::at for matrix of dimension %i by %i") % i % n % m;
        if ( j < 0 || j >= m ) throw formic::Exception("out-of-range col index j = %i in formic::Matrix::at for matrix of dimension %i by %i") % j % n % m;
        return this->_p[i+n*j];
      }

      const S & operator()(const L i, const L j) const {
        const int n = this->_n;
        const int m = this->_m;
        #ifndef NDEBUG
        if ( i < 0 || i >= n ) throw formic::Exception("out-of-range row index i = %i in formic::Matrix::operator() for matrix of dimension %i by %i") % i % n % m;
        if ( j < 0 || j >= m ) throw formic::Exception("out-of-range col index j = %i in formic::Matrix::operator() for matrix of dimension %i by %i") % j % n % m;
        #endif
        return this->_p[i+n*j];
      }

      S & operator()(const L i, const L j) {
        const int n = this->_n;
        const int m = this->_m;
        #ifndef NDEBUG
        if ( i < 0 || i >= n ) throw formic::Exception("out-of-range row index i = %i in formic::Matrix::operator() for matrix of dimension %i by %i") % i % n % m;
        if ( j < 0 || j >= m ) throw formic::Exception("out-of-range col index j = %i in formic::Matrix::operator() for matrix of dimension %i by %i") % j % n % m;
        #endif
        return this->_p[i+n*j];
      }

      S * begin() { return this->_p; }
      S * end()   { return this->_p + this->size(); }

      formic::Matrix<S,L> & operator=(const S a) {
        std::fill(this->_p, this->_p + this->_n * this->_m, a);
        return *this;
      }

      formic::Matrix<S,L> & operator+=(const S a) {
        for (L i = 0; i < this->size(); i++)
          this->_p[i] += a;
        return *this;
      }

      formic::Matrix<S,L> & operator-=(const S a) {
        for (L i = 0; i < this->size(); i++)
          this->_p[i] -= a;
        return *this;
      }

      formic::Matrix<S,L> & operator*=(const S a) {
        for (L i = 0; i < this->size(); i++)
          this->_p[i] *= a;
        return *this;
      }

      formic::Matrix<S,L> & operator/=(const S a) {
        for (L i = 0; i < this->size(); i++)
          this->_p[i] /= a;
        return *this;
      }

      formic::Matrix<S,L> & operator+=(const formic::ConstMatrix<S,L> & other) {
        if ( this->_n != other.rows() || this->_m != other.cols() )
          throw formic::Exception("cannot add matrices of different dimensions");
        for (L j = 0; j < this->_m; j++)
        for (L i = 0; i < this->_n; i++)
          (*this)(i,j) += other(i,j);
        return *this;
      }

      formic::Matrix<S,L> & operator-=(const formic::ConstMatrix<S,L> & other) {
        if ( this->_n != other.rows() || this->_m != other.cols() )
          throw formic::Exception("cannot subtract matrices of different dimensions");
        for (L j = 0; j < this->_m; j++)
        for (L i = 0; i < this->_n; i++)
          (*this)(i,j) -= other(i,j);
        return *this;
      }

      formic::Matrix<S,L> & axpy(const S a, const formic::ConstMatrix<S,L> & other) {
        if ( this->_n != other.rows() || this->_m != other.cols() )
          throw formic::Exception("cannot axpy matrices of different dimensions");
        for (L j = 0; j < this->_m; j++)
        for (L i = 0; i < this->_n; i++)
          (*this)(i,j) += a * other(i,j);
        return *this;
      }

      /// \brief  transpose the matrix in place
      formic::Matrix<S,L> & tip() {
        if (this->_n == this->_m) {
          for (L i = 0; i < this->_n; i++)
          for (L j = 0; j < i; j++)
            std::swap((*this)(i,j), (*this)(j,i));
        } else {
          const formic::ConstMatrix<S,L> temp = this->clone();
          std::swap(this->_n,this->_m);
          for (L j = 0; j < this->_m; j++)
          for (L i = 0; i < this->_n; i++)
            (*this)(i,j) = temp(j,i);
        }
        return *this;
      }

  };

}

template<class S, class L>
formic::Matrix<S,L> formic::ConstMatrix<S,L>::clone() const {
  formic::Matrix<S,L> retval(_n, _m);
  for (L j = 0; j < _m; j++)
  for (L i = 0; i < _n; i++)
    retval(i,j) = (*this)(i,j);
  return retval;
}

template<class S, class L>
formic::Matrix<S,L> formic::ConstMatrix<S,L>::operator-() const {
  formic::Matrix<S,L> retval(_n,_m);
  for (L j = 0; j < _m; j++)
  for (L i = 0; i < _n; i++)
    retval(i,j) = -(*this)(i,j);
  return retval;
}

template<class S, class L>
formic::Matrix<S,L> formic::ConstMatrix<S,L>::operator+() const {
  return this->clone();
}

/// \brief  return the matrix transpose in a new matrix
template<class S, class L>
formic::Matrix<S,L> formic::ConstMatrix<S,L>::t() const {
  formic::Matrix<S,L> retval(_m, _n);
  for (L i = 0; i < _n; i++)
  for (L j = 0; j < _m; j++)
    retval(j,i) = (*this)(i,j);
  return retval;
}

/// \brief  return the matrix conjugate transpose in a new matrix
template<class S, class L>
formic::Matrix<S,L> formic::ConstMatrix<S,L>::c() const {
  formic::Matrix<S,L> retval(_m, _n);
  for (L i = 0; i < _n; i++)
  for (L j = 0; j < _m; j++)
    retval(j,i) = formic::conj((*this)(i,j));
  return retval;
}

/// \brief  print the matrix
template<class S, class L>
std::string formic::ConstMatrix<S,L>::print(const std::string & fmt, const std::string & name) const {
  std::string retval;
  if ( name.size() > 0 )
    retval += "printing matrix \"" + name + "\":\n";
  for (L i = 0; i < _n; i++) {
    for (L j = 0; j < _m; j++)
      retval += ( boost::format("  " + fmt) % (*this)(i,j) ).str();
    retval += "\n";
  }
  return retval;
}

template<class S, class L>
template<class V>
formic::Matrix<S,L> formic::ConstMatrix<S,L>::slice(const V & row_indices, const V & col_indices) const {
  const int n = row_indices.size();
  const int m = col_indices.size();
  if ( n > _n )
    throw formic::Exception("too many row indices (%i) in formic::ConstMatrix::slice for matrix of dimension %i by %i") % n % _n % _m;
  if ( m > _m )
    throw formic::Exception("too many col indices (%i) in formic::ConstMatrix::slice for matrix of dimension %i by %i") % m % _n % _m;
  formic::Matrix<S,L> retval(n, m);
  for (L j = 0; j < m; j++)
  for (L i = 0; i < n; i++)
    retval(i,j) = (*this)(row_indices[i], col_indices[j]);
  return retval;
}

template<class S, class L> inline
formic::Matrix<S,L> operator+(const formic::ConstMatrix<S,L> & m1, const formic::ConstMatrix<S,L> & m2) {
  formic::Matrix<S,L> retval = m1.clone();
  retval += m2;
  return retval;
}

template<class S, class L> inline
formic::Matrix<S,L> operator-(const formic::ConstMatrix<S,L> & m1, const formic::ConstMatrix<S,L> & m2) {
  formic::Matrix<S,L> retval = m1.clone();
  retval -= m2;
  return retval;
}

template<class S, class L> inline
formic::Matrix<S,L> operator+(const formic::ConstMatrix<S,L> & m1, const S a) {
  formic::Matrix<S,L> retval = m1.clone();
  retval += a;
  return retval;
}

template<class S, class L> inline
formic::Matrix<S,L> operator-(const formic::ConstMatrix<S,L> & m1, const S a) {
  formic::Matrix<S,L> retval = m1.clone();
  retval -= a;
  return retval;
}

template<class S, class L> inline
formic::Matrix<S,L> operator*(const formic::ConstMatrix<S,L> & m1, const S a) {
  formic::Matrix<S,L> retval = m1.clone();
  retval *= a;
  return retval;
}

template<class S, class L> inline
formic::Matrix<S,L> operator/(const formic::ConstMatrix<S,L> & m1, const S a) {
  formic::Matrix<S,L> retval = m1.clone();
  retval /= a;
  return retval;
}



template<class S, class L> inline
formic::Matrix<S,L> operator+(const S a, const formic::ConstMatrix<S,L> & m1) {
  formic::Matrix<S,L> retval = m1.clone();
  retval += a;
  return retval;
}

template<class S, class L> inline
formic::Matrix<S,L> operator-(const S a, const formic::ConstMatrix<S,L> & m1) {
  formic::Matrix<S,L> retval = m1.clone();
  retval -= a;
  return retval;
}

template<class S, class L> inline
formic::Matrix<S,L> operator*(const S a, const formic::ConstMatrix<S,L> & m1) {
  formic::Matrix<S,L> retval = m1.clone();
  retval *= a;
  return retval;
}

template<class S, class L> inline
formic::Matrix<S,L> operator/(const S a, const formic::ConstMatrix<S,L> & m1) {
  formic::Matrix<S,L> retval(m1.rows(), m1.cols());
  for (L j = 0; j < m1.cols(); j++)
  for (L i = 0; i < m1.rows(); i++)
    retval(i,j) = a / m1(i,j);
  return retval;
}



template<class S, class L> inline
formic::Matrix<S,L> operator*(const formic::ConstMatrix<S,L> & m1, const formic::ConstMatrix<S,L> & m2) {
  if ( m1.cols() != m2.rows() )
    throw formic::Exception("cannot multiply: inner dimensions of (%i by %i) and (%i by %i) matrices do not match.")
          % m1.rows() % m1.cols() % m2.rows() % m2.cols();
  formic::Matrix<S,L> retval(m1.rows(), m2.cols(), formic::zero(S()));
  formic::xgemm('N', 'N', m1.rows(), m2.cols(), m1.cols(), formic::unity(S()),
                &m1.at(0,0), m1.rows(), &m2.at(0,0), m2.rows(), formic::zero(S()), &retval.at(0,0), retval.rows());
  return retval;
}

template<class S, class L> inline
formic::Vector<S,L> operator*(const formic::ConstMatrix<S,L> & m1, const formic::Vector<S,L> & v) {
  if ( m1.cols() != v.size() )
    throw formic::Exception("cannot multiply: dimensions of (%i by %i) matrix and (%i by 1) vector do not match.")
          % m1.rows() % m1.cols() % v.size();
  if ( v.size() == 0 )
    throw formic::Exception("cannot multiply: vector is zero length");
  if ( m1.rows() == 0 )
    throw formic::Exception("cannot multiply: matrix has no rows");
  return formic::Vector<S,L>(m1.rows()) <<= &( m1 * formic::ConstMatrix<S,L>(v.size(), 1, &v.at(0)) ).at(0,0);
}

template<class S, class L> inline
formic::Matrix<S,L> operator*(const formic::Vector<S,L> & v, const formic::ConstMatrix<S,L> & m1) {
  if ( m1.rows() != v.size() )
    throw formic::Exception("cannot multiply: dimensions of (1 by %i) vector and (%i by %i) matrix do not match.")
          % v.size() % m1.rows() % m1.cols();
  if ( v.size() == 0 )
    throw formic::Exception("cannot multiply: vector is zero length");
  if ( m1.cols() == 0 )
    throw formic::Exception("cannot multiply: matrix has no columns");
  return formic::Vector<S,L>(m1.cols()) <<= &( formic::ConstMatrix<S,L>(1, v.size(), &v.at(0)) * m1 ).at(0,0);
}

template<class S, class L>
formic::Matrix<S,L> formic::Vector<S,L>::t() const { return formic::Matrix<S,L>(_n, 1, _p).t(); }

template<class S, class L>
formic::Matrix<S,L> formic::Vector<S,L>::m() const { return formic::Matrix<S,L>(1, _n, _p).t(); }

namespace formic {

  template<class S, class L> inline
  formic::Matrix<S,L> identity_matrix(const L n) {
    formic::Matrix<S,L> retval(n, n, formic::zero(S()));
    for (L i = 0; i < n; i++)
      retval(i,i) = formic::unity(S());
    return retval;
  }

  template<class S, class L>
  formic::Matrix<S,L> matrix_exponent(const formic::ConstMatrix<S,L> & m,
                                      const double tol);

  template<class S, class L>
  formic::Matrix<S,L> matrix_exponent_der_adj(const formic::ConstMatrix<S,L> & m,
                                              const formic::ConstMatrix<S,L> & a,
                                              const double tol);

}

#endif
